import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import folium
import networkx as nx
import osmnx as ox
from pyproj import Proj, transform
from shapely.ops import nearest_points
from geopy.distance import geodesic
import pickle
import networkx as nx
from shapely.geometry import LineString
from scipy import spatial
from tqdm import tqdm
import gurobipy as gp
Location by Frequency¶
In this notebook, the location of the inspector is decided only by looking at the frequency of the accident in the road segment. The network shapefile is already provided road section id, hence it is going to be used to discritized the road.
The main work-flow for this method is:
- Calculate the distance between accident and the road network.
- Assign the accident according its road section id on the shapefiles.
- Assign the the road inspector to the highest accident frequency.
- Evaluate the travel time to each incident.
The calculation of distance between accident and road network¶
The first section is to calculate the distance. In this case, several data filtering is needed. The train and test data split are also done here. Other than that, the coordinate reference systems of the incident and the road network is synchronozed to the WGS84 by using the given algorithm and also existing python module (Pyproj)
# Functions from the teachers to transform the crs.
def DutchRDtoWGS84(rdX, rdY):
""" Convert DutchRD to WGS84
"""
RD_MINIMUM_X = 11000
RD_MAXIMUM_X = 280000
RD_MINIMUM_Y = 300000
RD_MAXIMUM_Y = 630000
if (rdX < RD_MINIMUM_X or rdX > RD_MAXIMUM_X
or rdY < RD_MINIMUM_Y or rdY > RD_MAXIMUM_Y):
resultNorth = -1
resultEast = -1
return resultNorth, resultEast
# else
dX = (rdX - 155000.0) / 100000.0
dY = (rdY - 463000.0) / 100000.0
k = [[3600 * 52.15517440, 3235.65389, -0.24750, -0.06550, 0.0],
[-0.00738 , -0.00012, 0.0 , 0.0 , 0.0],
[-32.58297 , -0.84978, -0.01709, -0.00039, 0.0],
[0.0 , 0.0 , 0.0 , 0.0 , 0.0],
[0.00530 , 0.00033, 0.0 , 0.0 , 0.0],
[0.0 , 0.0 , 0.0 , 0.0 , 0.0]]
l = [[3600 * 5.38720621, 0.01199, 0.00022, 0.0 , 0.0],
[5260.52916 , 105.94684, 2.45656, 0.05594, 0.00128],
[-0.00022 , 0.0 , 0.0 , 0.0 , 0.0],
[-0.81885 , -0.05607, -0.00256, 0.0 , 0.0],
[0.0 , 0.0 , 0.0 , 0.0 , 0.0],
[0.00026 , 0.0 , 0.0 , 0.0 , 0.0]]
resultNorth = 0
resultEast = 0
powX = 1
for p in range(6):
powY = 1
for q in range(5):
resultNorth = resultNorth + k[p][q] * powX * powY / 3600.0
resultEast = resultEast + l[p][q] * powX * powY / 3600.0
powY = powY * dY
powX = powX * dX
return resultNorth, resultEast
def WGS84toDutchRD(wgs84East, wgs84North):
# translated from Peter Knoppers's code
# wgs84East: longtitude
# wgs84North: latitude
# Western boundary of the Dutch RD system. */
WGS84_WEST_LIMIT = 3.2
# Eastern boundary of the Dutch RD system. */
WGS84_EAST_LIMIT = 7.3
# Northern boundary of the Dutch RD system. */
WGS84_SOUTH_LIMIT = 50.6
# Southern boundary of the Dutch RD system. */
WGS84_NORTH_LIMIT = 53.7
if (wgs84North > WGS84_NORTH_LIMIT) or \
(wgs84North < WGS84_SOUTH_LIMIT) or \
(wgs84East < WGS84_WEST_LIMIT) or \
(wgs84East > WGS84_EAST_LIMIT):
resultX = -1
resultY = -1
else:
r = [[155000.00, 190094.945, -0.008, -32.391, 0.0],
[-0.705, -11832.228, 0.0 , 0.608, 0.0],
[0.0 , -114.221, 0.0 , 0.148, 0.0],
[0.0 , -2.340, 0.0 , 0.0 , 0.0],
[0.0 , 0.0 , 0.0 , 0.0 , 0.0]]
s = [[463000.00 , 0.433, 3638.893, 0.0 , 0.092],
[309056.544, -0.032, -157.984, 0.0 , -0.054],
[73.077, 0.0 , -6.439, 0.0 , 0.0],
[59.788, 0.0 , 0.0 , 0.0 , 0.0],
[0.0 , 0.0 , 0.0 , 0.0 , 0.0]]
resultX = 0
resultY = 0
powNorth = 1
dNorth = 0.36 * (wgs84North - 52.15517440)
dEast = 0.36 * (wgs84East - 5.38720621)
for p in range(5):
powEast = 1
for q in range(5):
resultX = resultX + r[p][q] * powEast * powNorth
resultY = resultY + s[p][q] * powEast * powNorth
powEast = powEast * dEast
powNorth = powNorth * dNorth
return resultX, resultY
def calc_distance(line_wkt):
line = ogr.CreateGeometryFromWkt(line_wkt)
points = line.GetPoints()
d = 0
for p0, p1 in zip(points, points[1:]):
d = d + geodesic(p0, p1).m
return d
if __name__=="__main__":
x, y = WGS84toDutchRD(4.33, 52.04)
print(DutchRDtoWGS84(x, y))
(52.03999999894767, 4.330000046074026)
# Function from the pyproj module to project the network shapefiles entirely
def utm_projection(gdf,initial_crs):
target_crs = 'EPSG:32631'
transformer = Proj(init=initial_crs), Proj(init=target_crs)
gdf_utm = gdf.to_crs(target_crs)
return(gdf_utm)
def wgs84_projection(gdf, initial_crs):
target_crs = 'EPSG:4326'
transformer = Proj(init=initial_crs), Proj(init=target_crs)
gdf_wgs84 = gdf.to_crs(target_crs)
return(gdf_wgs84)
In the cell below, the network data shapefile is loaded and transformed.
# Extract subnetwork
highway_shapefile = 'Shapefiles/Snelheid_Wegvakken.shp'
network_temp = gpd.read_file(highway_shapefile)
# Make a geodataframe of the network and
# define the crs to make the transformation easier
df_rd = gpd.GeoDataFrame(geometry=network_temp['geometry'], crs='EPSG:28992')
df_utm = utm_projection(df_rd, 'EPSG:28992')
df_wgs84 = wgs84_projection(df_rd, 'EPSG:28992')
# Changing the column to english for easier analysis
df_utm['Road_section_id'] = network_temp['WVK_ID']
df_utm['Road_number'] = network_temp['WEGNR_HMP']
df_wgs84['Road_section_id'] = network_temp['WVK_ID']
df_wgs84['Road_number'] = network_temp['WEGNR_HMP']
c:\Users\syaka\anaconda3\envs\geospatial\Lib\site-packages\pyproj\crs\crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string) c:\Users\syaka\anaconda3\envs\geospatial\Lib\site-packages\pyproj\crs\crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string) c:\Users\syaka\anaconda3\envs\geospatial\Lib\site-packages\pyproj\crs\crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string) c:\Users\syaka\anaconda3\envs\geospatial\Lib\site-packages\pyproj\crs\crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string)
The incident data frame is loaded and filtered for only the highway incidents. In this part, the train-test data split is also made.
# Load incident csv data
incident = pd.read_csv('incidents19Q3Q4.csv')
# Rename the column for easier analysis
incident.rename(columns = {'vild_primair_wegnummer':'Road_number',
'primaire_locatie_lengtegraad':'longitude',
'primaire_locatie_breedtegraad':'latitude'}, inplace=True)
# Seperate the highway data
incident = incident.loc[np.where(incident['Road_number'].str[0] == 'A')]
# Train-test data split
data = incident[['longitude', 'latitude']].values
num_samples = len(data)
num_samples_keep = int(0.8*num_samples)
np.random.shuffle(data)
train_data = incident.iloc[:num_samples_keep]
test_data = incident.iloc[num_samples_keep:]
train_data
| Unnamed: 0 | id | type | starttime_new | endtime_new | Road_number | longitude | latitude | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0 | A-ALL-IM19087676NLD_1 | vehicle_obstruction | 2019-08-28 12:11:32 | 2019-12-11 11:32:28 | A1 | 4.974663 | 52.346931 |
| 1 | 3 | A-ALL-IM19087677NLD_1 | vehicle_obstruction | 2019-08-28 12:11:32 | 2019-12-11 11:32:28 | A9 | 4.716725 | 52.514820 |
| 2 | 5 | A-ALL-IM19087678NLD_1 | vehicle_obstruction | 2019-08-28 12:11:32 | 2019-12-11 11:32:28 | A9 | 4.738364 | 52.609730 |
| 3 | 7 | A-ALL-IM19087679NLD_1 | vehicle_obstruction | 2019-08-28 12:11:32 | 2019-12-11 11:32:28 | A35 | 6.824692 | 52.204929 |
| 4 | 174938 | A-ALL-IM19087680NLD_1 | vehicle_obstruction | 2019-08-28 12:11:32 | 2019-12-11 11:32:28 | A4 | 4.346407 | 52.041920 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 71333 | 337797 | RWS03_804572_1 | vehicle_obstruction | 2019-10-29 22:19:50 | 2019-10-29 23:35:46 | A2 | 5.791368 | 51.003658 |
| 71334 | 337800 | RWS03_804573_1 | general_obstruction | 2019-10-29 22:32:55 | 2019-10-29 23:02:59 | A15 | 4.372907 | 51.875309 |
| 71335 | 337803 | RWS03_804574_1 | accident | 2019-10-29 22:41:16 | 2019-10-29 23:56:46 | A1 | 6.343580 | 52.244808 |
| 71336 | 337804 | RWS03_804575_1 | accident | 2019-10-29 22:49:03 | 2019-10-29 23:12:35 | A9 | 4.705645 | 52.396221 |
| 71337 | 337806 | RWS03_804577_1 | vehicle_obstruction | 2019-10-29 23:03:09 | 2019-10-30 00:18:46 | A13 | 4.351567 | 52.038780 |
57690 rows × 8 columns
Based on the train data, the geodataframe of the incident is made.
# Incident geodataframe from the train dataset
Geometry = gpd.points_from_xy(train_data.latitude, train_data.longitude)
gdf_incident_wgs84 = gpd.GeoDataFrame(train_data, geometry=Geometry, crs='EPSG:4326')
gdf_incident_utm = utm_projection(gdf_incident_wgs84, 'EPSG:4326')
c:\Users\syaka\anaconda3\envs\geospatial\Lib\site-packages\pyproj\crs\crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string) c:\Users\syaka\anaconda3\envs\geospatial\Lib\site-packages\pyproj\crs\crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string)
# Incident heatmap visualisation
from folium import plugins
map = folium.Map(location=[52.399190, 4.893658])
gjson = folium.features.GeoJson(
df_wgs84,
).add_to(map)
# List of accident for the heatmap visualisation
geo_df_list = [[gdf_incident_wgs84['geometry'].iloc[i].xy[0][0],
gdf_incident_wgs84['geometry'].iloc[i].xy[1][0]] for i in range(len(gdf_incident_wgs84.geometry))]
incident_a200_mark = plugins.HeatMap(geo_df_list).add_to(map)
map